#########################################################################################
# 
# Gender Specific Models
#
#########################################################################################

#########################################################################################
### Democratic Republic of Congo

drc.out.f <- ictreg.joint(Y ~ age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                            + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev,  
                            treat="treatment", 
                            outcome="soc_part.2",
                            outcome.reg="logistic",
                            constrained=TRUE,
                            J=3, data=congo[congo$female==1,]) # Model for Females
summary(drc.out.f)

drc.out.m <- ictreg.joint(Y ~ age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                            + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev,  
                            treat="treatment", 
                            outcome="soc_part.2",
                            outcome.reg="logistic",
                            constrained=TRUE,
                            J=3, data=congo[congo$female!=1,]) # Model for Males
summary(drc.out.m)

sim.eff.fun.gen <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$par.outcome,  m.out$vcov[25:37,25:37]) # Informal Bayes
  
  x.0 <- cbind(m.out$x, rep(0, dim(m.out$x)[1]))
  x.1 <- cbind(m.out$x, rep(1, dim(m.out$x)[1]))
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

drc.f <- round(mean(sim.eff.fun.gen(drc.out.f)), 2)
drc.ci.f <- round(quantile(sim.eff.fun.gen(drc.out.f), c(.05, .95)), 2)

drc.m <- round(mean(sim.eff.fun.gen(drc.out.m)), 2)
drc.ci.m <- round(quantile(sim.eff.fun.gen(drc.out.m), c(.05, .95)), 2)

drc.diff <- round(mean(sim.eff.fun.gen(drc.out.m)-sim.eff.fun.gen(drc.out.f)), 2)
drc.ci.diff <- round(quantile(sim.eff.fun.gen(drc.out.m)-sim.eff.fun.gen(drc.out.f), c(.05, .95)), 2)

drc.out.diff <- drc.out.f$par.outcome-drc.out.m$par.outcome
drc.out.diff.se <- sqrt(drc.out.f$se.outcome^2 + drc.out.m$se.outcome^2)

#########################################################################################
### Liberia
# ATTENTION: Bayesian simulation may take a while! 

liberia.f <- liberia[liberia$female==1,]
liberia.m <- liberia[liberia$female!=1,]

mod10.f = bayeslist(Y ~ age.z + edu.z + income.z + hh_size.z + cw_kill + county.1 + county.2, 
                    data = liberia.f,
                    treat = "treatment",
                    J = 3,
                    type = "predict",
                    predictvar = "outcome_ca",
                    predictvar_type = "binary",
                    nsim = 5000,
                    thin = 1,
                    CIsize = 0.95,
                    nchain = 2,
                    seeds = 342321,
                    prior = NULL,
                    parallel = T) # Model for Females

summary(mod10.f)

mod10.m = bayeslist(Y ~ age.z + edu.z + income.z + hh_size.z + cw_kill + county.1 + county.2, 
                    data = liberia.m,
                    treat = "treatment",
                    J = 3,
                    type = "predict",
                    predictvar = "outcome_ca",
                    predictvar_type = "binary",
                    nsim = 5000,
                    thin = 1,
                    CIsize = 0.95,
                    nchain = 2,
                    seeds = 342321,
                    prior = NULL,
                    parallel = T) # Model for Females

summary(mod10.m)

#lib.out.f <- readRDS("Submission/AJPS/AJPS_replication_materials/lib_out_f.rds")
#lib.out.m <- readRDS("Submission/AJPS/AJPS_replication_materials/lib_out_m.rds")

# Effect Size Female
x.0.f <- cbind(lib.out.f$X, rep(0, dim(lib.out.f$X)[1]))
x.1.f <- cbind(lib.out.f$X, rep(1, dim(lib.out.f$X)[1]))

pred.0.f <- invlogit(x.0.f%*%t(lib.out.f$sampledf[,18:26]))
pred.1.f <- invlogit(x.1.f%*%t(lib.out.f$sampledf[,18:26]))

pred.eff.f <- pred.1.f - pred.0.f
pred.eff.av.f <- apply(pred.eff.f, 2, mean)

lib.f <- round(mean(pred.eff.av.f), 2)
lib.ci.f <- round(quantile(pred.eff.av.f, c(.05, .95)), 2)

# Effect Size Male
x.0.m <- cbind(lib.out.m$X, rep(0, dim(lib.out.m$X)[1]))
x.1.m <- cbind(lib.out.m$X, rep(1, dim(lib.out.m$X)[1]))

pred.0.m <- invlogit(x.0.m%*%t(lib.out.m$sampledf[,18:26]))
pred.1.m <- invlogit(x.1.m%*%t(lib.out.m$sampledf[,18:26]))

pred.eff.m <- pred.1.m - pred.0.m
pred.eff.av.m <- apply(pred.eff.m, 2, mean)

lib.m <- round(mean(pred.eff.av.m), 2)
lib.ci.m <- round(quantile(pred.eff.av.m, c(.05, .95)), 2)

lib.diff <-  round(mean(pred.eff.av.m-pred.eff.av.f), 2)
lib.ci.diff <-round(quantile(pred.eff.av.m-pred.eff.av.f, c(.05, .95)), 2)


#########################################################################################
### Sri Lanka

sri.out.f <- ictreg.joint(Y ~  age.z + edu.z + income.z + hh_size.z +  killed + trauma +  prov.8 + prov.9 + prior,
                            treat="treatment", 
                            outcome="soc_part.2",
                            outcome.reg="logistic",
                            constrained=TRUE,
                            J=3, data=sri[sri$female==1,])
summary(sri.out.f)

sri.out.m <- ictreg.joint(Y ~  age.z + edu.z + income.z + hh_size.z +  killed + trauma +  prov.8 + prov.9 + prior ,
                            treat="treatment", 
                            outcome="soc_part.2",
                            outcome.reg="logistic",
                            constrained=TRUE,
                            J=3, data=sri[sri$female!=1,])
summary(sri.out.m)

sim.eff.fun.gen <- function(m.out){ 
  coef.sim <- mvrnorm(10000, m.out$par.outcome,  m.out$vcov[21:31,21:31]) # Informal Bayes
  
  x.0 <- cbind(m.out$x, rep(0, dim(m.out$x)[1]))
  x.1 <- cbind(m.out$x, rep(1, dim(m.out$x)[1]))
  
  pred.0 <- invlogit(x.0%*%t(coef.sim))
  pred.1 <- invlogit(x.1%*%t(coef.sim))
  
  pred.eff <- pred.1 - pred.0
  pred.eff.av <- apply(pred.eff, 2, mean)
  return(pred.eff.av)
}

sri.f <- round(mean(sim.eff.fun.gen(sri.out.f)), 2)
sri.ci.f <- round(quantile(sim.eff.fun.gen(sri.out.f), c(.05, .95)), 2)

sri.m <- round(mean(sim.eff.fun.gen(sri.out.m)), 2)
sri.ci.m <- round(quantile(sim.eff.fun.gen(sri.out.m), c(.05, .95)), 2)

sri.diff <- round(mean(sim.eff.fun.gen(sri.out.m)-sim.eff.fun.gen(sri.out.f)), 2)
sri.ci.diff <- round(quantile(sim.eff.fun.gen(sri.out.m)-sim.eff.fun.gen(sri.out.f), c(.05, .95)), 2)

m.out.1.a.diff <- sri.out.f$par.outcome-sri.out.m$par.outcome
m.out.1.a.diff.se <- sqrt(sri.out.f$se.outcome^2 + sri.out.m$se.outcome^2)

###############################################################################################

